MM = 100000;

defpb5

c = zeros(2*n+M,1);
A= zeros(3*n+2*M,2*n+M);
b = zeros(3*M+2*n,1);
Aeq = [];
beq = [];
lb = [-MM; 0];
ub = [MM; 1];
yidx = [false(2*n, 1); true(M, 1)];

% calcul de c
for i = n+1 : 2*n
	c(i) = 1;
end

% Calcul de A
for i = 1:n
	A(i,i)=1;
	A(i,i+n)=-1;
end
for i = (n+1) : (2*n)
	A(i, i-n) = -1;
	A(i, i-n+n) = -1;
end
for i = (2*n+1): 3*n
	for j = 1 :(i-2*n)
		A(i,j) = -(i-2*n)+j-1;
	end
end

for i = 3*n+1:3*n+M
    	A(i, i-n) = -MM;
end
for i = 3*n+M+1:3*n+2*M
    	A(i, i-n-M) = MM;
end

for i = 3*n+1:3*n+M
	timp = timpact(i-3*n);
	for j = 1:n
		if timp ~= 0 
			A(i,j) = timp;
			timp = timp - 1;
		else
			A(i,j) = 0;
		end
	end
end
for i = 3*n+1+M:3*n+M+M
	timp = timpact(i-3*n-M);
	for j = 1:n
		if timp ~= 0 
			A(i,j) = -timp;
			timp = timp - 1;
		else
			A(i,j) = 0;
		end
	end
end

% Calcul de b
for i=1:(2*n)
	b(i) = 0;
end
for i=(2*n+1):3*n
	b(i) = (i-2*n)*v0;
end
for i = (3*n+1): 3*n+M
	b(i) = pimmin(i - 3*n)- b(i-M);
end
for i = (3*n+M+1) : 3*n + 2*M
	b(i) = MM - pimmax(i -3*n-M)+b(i-2*M);
end

A

